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Abstract 



^ ' Solutions to Monge-Kantorovich equations, expressing optimality condition in 

■ mass transportation problem with cost equal to distance, are stationary points of a 
i critical-slope model for sand surface evolution. Using a dual variational formulation 
' of sand model, we compute both the optimal transport density and Kantorovich 
' potential as t — > oo limit of evolving sand flux and sand surface, respectively. 

o 

■ 1 Introduction 

■ 

S ' The Monge-Kantorovich equations, 

> 

/+-/- = -V-(aVn), (1) 
|Vm| < 1, a>0, |Vm| < 1 — >a = 0, (2) 

appear as a condition of optimality in the classical Monge-Kantorovich problem of optimal 
mass transportation with the cost function equal to the distance of transportation ^1 
121 El- Here are two given distributions of mass, a is the transport density, and u 
is the Kantorovich potential. Let A^(i?"') be the set of bounded Radon measures on 
i?", M+{R"') the subset of nonnegative measures, and Lip]^(i?"') the set of real valued 
Lipschitz functions with Lipschitz constant not greater than 1. If G A^+(-R") and 
satisfy f^{R^) = /~(-R"), there exist a transport density and a corresponding potential, 
{a, u} G M+{R'^) X Lip]^(_R"'), satisfying (P)-© in a weak sense, see p. If either /+ or f~ 
are absolutely continuous with respect to Lebesgue measure in the transport density 
a is unique and also absolutely continuous fTj. Obviously, the Kantorovich potential u is 
determined up to an additive constant. If the constant is fixed, the potential is still not 
unique outside spt(a). Note that the system (ll}-© arises also in a mass optimization 
problem 
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A critical-slope model for sand surface evolution ([S], see also [HI EI) 

f - dtu = -V ■ (aVu), u\t=o = uo{x), (3) 
|Vm| <k, a > 0, |Vm| <k — > a = 0, (4) 

can be regarded as a non-stationary version of the Monge-Kantorovich equations. Here 
X e R^, u{x,t) is the evolving sand surface, f{x,t) is the intensity of external source, k is 
the tangent of the sand angle of repose, a{x, t) is the sand transport density (if = 1) or 
proportional to it, and the initial surface Uq satisfies the equilibrium constraint | VmqI < k. 

Indeed, it can be proved that if the source intensity / is nonnegative, the term dtu 
is nonnegative too and, as was noted in [H:, comparison of the two systems shows that 
at each time moment sand from the source is instantly incorporated into the bulk in a 
way that minimizes sand transport under the constraint |Vu| < k. The height function 
u plays here the role of Kantorovich potential. 

In this note we explore an opposite approach to this similarity. Let the source intensity 
be a time-independent Radon measure, / = — f~, such that G A^+(i?") and satisfy 
f~^{R"') = f~{R^). Starting with an admissible initial condition uq, we seek a solution to 
Monge-Kantorovich equations as the t — > cxd limit of solution to Q-©- To solve this latter 
evolutionary problem we employ the numerical algorithm P based on a dual variational 
formulation of the sand model written for flux of sand pouring down the evolving sand 
surface [TU]. 



2 Variational formulations of sand model 

To study the evolution of sand surface, it is convenient to exclude from equations ©-Q 
the transport density a, the Lagrange multiplier related to the constraint |Vu| < k, and 
rewrite the model as a variational inequality for u alone. For bounded domains and ap- 
propriate boundary conditions, the inequality was obtained in jSlini- Similar formulation 
of the initial value problem has been independently derived and studied in [Zj: 

u{.,t)eK: (dtU- f,ip-u)>0, \/ipeK, 
u\t=o = Mo, 

where K = {ip E L?{R^) : |Vv?| < k a.e.}. We assume both Uq E K and / have compact 
supports. In this case spt(M) remains compact for each t > 0, see jTj. 

The inequality (0) can be used for efficient computation of the evolving sand surface 
u. However, the dual variable, a, remains unknown and difficult to find. Such situation 
is typical also of other critical-state problems. Because of this reason the dual variational 
formulations resembling mixed variational inequalities in elasto-plasticity JT] and 
allowing to find both variables simultaneously (see [H]), can be preferable. Having in 
mind the application to Monge-Kantorovich problem, we assume x G -R" and allow for 
non-constant coefficients k = k{x) > to treat problems with, say, subregions through 
which mass transportation is forbidden or, on the contrary, where it is free of charge. 

Suppose the pair {u,a} satisfies the model relations ©-(HI)- Then, for the horizontal 
projection of surface fiux of sand q = —aVu and arbitrary test field '0 we obtain 

Vu ■ {ip - q) > -|Vn||V'| - Vu-q = -|Vu||V'| + k\q\ > -k\ip\ + k\q\. 
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Noting that support of u is bounded and integrating, we obtain (Vn, i/^ — q) > (j){q)—(j){il^), 
where 0(q) = /^^ k\q\, and also (Vu, il> — q) = — (m, V ■ {i/' — q})- Therefore, 

0(V)-0(g)-KV-{V-g})>O. 

Integrating now the balance equation ^ in time we get 

u = Uo + [ fdt-V ■ [ qdt (6) 
Jo Jo 

and arrive at the mixed variational inequality written for the sand flux alone: 

q{.,t)eV: (v-j^qdt-J^,V-{i^-q}^+<P{'4^)-<P{q)>0 (7) 

for any ij^ eV . Here JF = mq + f dt and we define 

V = e {M{R')X' : V- V G L2(i?")}. 

Several comments about the variational formulation ((7j) can be necessary. 

The problem is not coercive in the usual Sobolev spaces. Nevertheless, under suitable 
conditions on the source function /, existence of a solution can be established 

Since the divergence of flux q belongs to L^, we expect the continuity of normal com- 
ponent of sand flux; the tangential component of q can be discontinuous (this corresponds 
to continuity of transport density along the transport rays). To solve such a problem nu- 
merically, one should use the divergence-conforming finite elements that do not enforce 
any additional smoothness on the solution. 

Suppose the flux has been found as a solution to ((Zj). Obviously, the transport density 
is easily determined for A; > as a = \q\/k. The surface u (potential in the Monge- 
Kantorovich problem) can be calculated by means of the equation ©. 



3 Solution of Monge-Kantorovich equations 

To approximate the inequality (|7j) numerically, we smoothed the non-differentiable func- 
tional by introducing \q\e = {Iq]"^ + s'^Y^'^ , discretized the regularized equality problem 
in time, employed the divergence- conforming Raviart-Thomas finite elements of lowest 
order (see, e.g., ^3]) and used vertex sampling on the nonlinear term. Resulting nonlin- 
ear algebraic systems were solved at each time level iteratively using a form of successive 
over-relaxation. We refer to [12] for details and make only several remarks related to 
efficient implementation of this algorithm to Monge-Kantorovich equations. 

i) The needed solutions are obtained as the t ^ oo limit of solutions to ((Tj) with the 
time-independent sources satisfying /(-R") = 0. To find this limit, we need not solve 
the nonlinear equations at each time level with high accuracy: only a few iterations are 
needed before the transition to a new time level. If, however, the Kantorovich potential is 
also of interest by some reason, the time step and the number of iterations should ensure 
accurate calculation of the integral in (0) . 
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ii) Suppose that transportation through an open domain Q is forbidden. To model this 
situation, it is possible to solve the problem in R'^\fl with qn\dQ = 0. Another possibility 
(used in this work) is to solve also in Q but set k\Q = oo. Then, if q solves (|7j), the flux 
qIq must be zero. Since the normal component of flux is continuous, only tangential to 
domain boundary flux of mass is permitted: the optimal mass transportation rounds the 
obstacle. 

iii) Numerical solution of is possible only in a bounded domain and it is desirable 
to make the computational domain as small as possible. However, if the domain is not 
large enough, at some (depending on uq) time moment spt{q{.,t)) can reach the domain 
boundary. To avoid any material loss through the boundary, we can surround the do- 
main by an obstacle (see previous comment). On the other hand, we do not want the 
artificial obstacle to affect our solution. Since the transport set in Monge-Kantorovich 
problem (without obstacles) consists of straight transport rays that begin in spt(/"*") and 
end in spt(/~), we can choose any computational domain containing the convex hull of 
spt(/) =spt(/"*") IJ spt(/~) and surround it by an obstacle. Our choice of domain (and 
Uq) can influence only the Kantorovich potential outside the establishing in t — > oo limit 
transport set, where this potential is not unique. 

In all examples below, the computations were performed in a square surrounded by an 
obstacle (we set k = 10^ in a thin border around the square). Matlab Partial Differential 
Equation Toolbox was used for generation of finite element meshes. 
Example 1. Support of /"*" consists of two upper ellipses in which the density of is 1 
(see Fig. [TJ. The mass from /"•" is to be uniformly distributed in the third ellipse below. 




Figure 1: Solution to Monge-Kantorovich equations. Left: the optimal flux and levels of 
Kantorovich potential. Right: contour plot of \q\. 

Example 2. Let and /~ be uniformly distributed in upper and lower rectangles, 
respectively, the ellipse between them is an obstacle (see Fig. E)). In this case the transport 
density is not absolutely continuous: it becomes concentrated on part of the obstacle 
boundary (see also for analysis of transport density in a sandpile growth problem 
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with an obstacle). 




Figure 2: As in fig. d problem with an obstacle to mass transportation. 

Example 3. Let and /~ be uniformly distributed in upper and lower rectangles, 
respectively, and the third polygon in Fig. |21 be a "highway" where the transportation 
cost is significantly reduced. To model this situation, we set k = 0.01 in the polygon 
(which makes transportation through this area hundred times cheaper) . Not surprisingly, 
the optimal transportation in this case is concentrated mostly upon the highway. 




Figure 3: As in fig. ^ problem with a region of cheap transportation. 



4 Conclusion 



We used formally the similarity of the Monge-Kantorovich equations and a model for 
sand surface evolution to present solutions to Monge-Kantorovich equations as stationary 
points of a dynamical system and find them numerically. We note that, although it 
remains to prove rigorously the solutions to non-stationary problem do have the t —>■ oo 
limit, numerical examples confirm this statement and demonstrate the efficiency of our 
approach. Our hypothesis is the convergence, in fact, occurs in a finite time. 

It seems interesting to give also an "optimal mass transportation" interpretation of 
sand model equations. Let f^{x,t) and f~{x,t) be, respectively, the production and 
consumptions rates of some homogeneous commodity and uo{x) its initial distribution. 
Let, at any time moment, the local price of this commodity be defined as —u{x,t), i.e., 
the price is positive if there is a deficit, u{x,t) < 0, and is negative if there is a surplus 
u{x,t) > which needs storing. 

Whenever the price difference between two arbitrary points, x and y, exceeds the cost 
of transportation of one commodity unit from one of these points to another, c{x, y) = 
\x — y\, the prices adjust themselves instantaneously, since a cheaper price is available by 
buying elsewhere and transporting. Condition of price equilibrium can thus be written in 
the form \u{x) — u{y)\ < \x — y\ or, locally, |Vu| < 1. Assuming the initial distribution 
of goods uo satisfies this condition, we describe its further evolution driven by production 
and consumption as 

dtu + V ■ q = - u\t=o = Uo{x) 

and assume that the mass fiux q is always directed towards the price gradient: q = —aVu, 
where a{x,t) > is an unknown scalar function. Finally, if \u{x) — u{y)\ < |a; — y\, the 
transportation from one of these points to another can only increase the cost and is 
unprofitable. Locally, this can be reformulated as the condition |Vm| < 1 — > a{x,t) = 0, 
which brings us to the critical-state model ©-© in which all transport occurs at the 
border of equilibrium. 
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